Particle Set

Simulate trajectories of a particle cloud in a two-dimensional flow field. A doubly-periodic domain and randomly-generated flow fields are initially used. For additional documentation e.g. see : 1, 2, 3, 4

Exercises:

• change the initial distribution of partices
• increase the duration of the trajectories simulation
• treat the non-periodic domain case by padding u,v with zeros
• replace u,v with your own two-dimensional flow fields

1. Import Software

using IndividualDisplacements, DataFrames
p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"../examples/flow_fields.jl"));

2. Flow Fields

A convenient way to set up the flow fields using the MeshArrays.jl package (which handles such staggered grids in general fashion) is to call the convert_to_FlowFields function with u,v arrays and time span as arguments.

u,v,ϕ=random_flow_field()

𝐹=convert_to_FlowFields(u,v,10.0);

In the above, the u,v arrays generated by random_flow_field can be replaced with any other pair provided by the user.

A couple of important considerations, however:

• u,v are staggered on a C-grid; by -0.5 grid point in direction 1 for u (2 for v)

from the grid cell center (0.5,0.5)

• u,v here derive from streamfunction ϕ, defined at the corner point, which ensures that

the resulting u,v is non-divergent, purely rotational, over the C-grid domain. In brief:

u=-(circshift(ϕ, (0,-1))-ϕ)
v=(circshift(ϕ, (-1,0))-ϕ)

If user were to start with collocated velocity (uC,vC at the grid cell center) then one can easily obtain the staggered velocity (u,v) as follows. These may contain both rotational and divergent components.

u=0.5*(circshift(uC, (0,1))+uC)
v=0.5*(circshift(vC, (1,0))+vC)

3. Initialize Individuals

For example, we can initialize 100 particles within a central subdomain as follows.

np,nq=size(u)
x=np*(0. .+ 1.0*rand(1000))
y=nq*(0. .+ 1.0*rand(1000))
a=ones(size(x)); #subdomain array index (just 1 here)

The following constructor function wraps everything in the Individuals data structure.

𝐼=Individuals(𝐹,x,y,a)

  📌 details     = (1, 1000) Array{Float64,1}
  🔴 details     = (0, 5) ["ID", "x", "y", "fid", "t"]
  🆔 range       = (1, 1000)
  🚄 function    = dxy_dt!
  ∫  function    = default_solver
  🔧 function    = postprocess_MeshArray
  𝑃  details     = (:u0, :u1, :v0, :v1, :𝑇, :update_location!)

3. Compute Trajectories

The time period is 𝐼.𝑃.𝑇 by default, unless ∫!(𝐼,𝑇) is called instead.

∫!(𝐼)

1×1000 Array{Array{Float64,1},2}:
 [7.27738, 13.0045, 1.0]  [4.47776, 11.5454, 1.0]  …  [5.01068, 6.29591, 1.0]

4. Plot Results

For example, generate a simple animation:

p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"../examples/recipes_plots.jl"));

🔴_by_t = groupby(𝐼.🔴, :t)
anim = @animate for t in eachindex(🔴_by_t)
   phi_scatter(ϕ,🔴_by_t[t])
end

pth=tempdir()*"/"
gif(anim, pth*"RandomFlow.gif", fps = 10)

This page was generated using Literate.jl.